%% "Econometrics of the Hodrick-Prescott Filter," Robert de Jong and Neslihan Sakarya (2015)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% This program compiles Figure 1 (in Appendix 1) of the paper

%%

clear all;
close all;
clc;
T=50; %% total number of observations
lm=1600; %% the smoothing parameter, lambda


%% Defining function f_T_lambda, it is denoted as f

a=zeros(T,T); % the matrix for cos(pi*(j-1)*m/T) // used in f function 
df=zeros(T,1); % the vector for 1/(1+16*lm*sin(pi*(j-1)/(2*T))^4) // used in both f 

for j=2:T;
    for m=1:2*T;
        a(m,j)=cos(pi*(j-1)*m/T);
        df(j)=1/(1+16*lm*sin(pi*(j-1)/(2*T))^4);
    end
end
b=a*df; %% sum j=2 to T cos(pi*(j-1)*m/T)*(1+16*lm*sin(pi*(j-1)/(2*T))^4)^(-1)

f=zeros(T,1); 

for m=1:T; 
   f(m)=(2*T)^(-1)+cos(pi*m)*(2*T)^(-1)*(1+16*lm)^(-1)+T^(-1)*b(m); %% f function and note that cos(pi*m)=(-1)^m
end

%% This part is for defining f(0) since there is no zero entry of the vector f, it is called f0 and all the corresponding variables will be defined similarly, 
 % and note that cos(pi*(j-1)*m/T)=1 when m=0

 b0=0;
for j=2:T;
    b0=b0+(1+16*lm*sin(pi*(j-1)/(2*T))^4)^(-1); % sum j=2 to T (1+16*lm*sin(pi*(j-1)/(2*T))^4)^(-1)
end

f0=(2*T)^(-1)+(2*T)^(-1)*(1+16*lm)^(-1)+T^(-1)*b0; % Note that cos(pi*m)=1 when m=0

%% The results below will be displayed on the command window

ff=[f0; f] %this vector gives the value of f_T_lambda for m=0,1,2,...,T
printmat(ff, 'f_T_lambda'); %the values of f function with the name f_T_lambda

m=-T:1:T; 
f_T_lambda=[flipud(f); f0; f]; %% note that flipud(f)=f since f_T_lambda is a symmetric function in other words, f_T_lambda(m)=f_T_lambda(-m)
plot(m, f_T_lambda)
title('The Function f')
xlabel('m')
ylabel('Function Values')